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Abstract 

Sexually transmitted infections (STIs) continue to present a complex and costly challenge to public 
health programs. The preferences and social dynamics of a population can have a large impact on 
the course of an outbreak as well as the effectiveness of interventions intended to influence individ¬ 
ual behavior. In addition, individuals may alter their sexual behavior in response to the presence of 
STIs, creating a feedback loop between transmission and behavior. We investigate the consequences 
of modeling the interaction between STI transmission and prophylactic use with a model that links 
a Susceptible-Infectious-Susceptible (SIS) system to evolutionary game dynamics that determine the 
effective contact rate. The combined model framework allows us to address protective behavior by both 
infected and susceptible individuals. Feedback between behavioral adaptation and prevalence creates a 
wide range of dynamic behaviors in the combined model, including damped and sustained oscillations 
as well as bistability, depending on the behavioral parameters and disease growth rate. We found that 
disease extinction is possible for multiple regions where TZo > 1 , due to behavior adaptation driving 
the epidemic downward, although conversely endemic prevalence for arbitrarily low TZo is also possible 
if contact rates are sufficiently high. We also tested how model misspecification might affect disease 
forecasting and estimation of the model parameters and TZq. We found that alternative models that 
neglect the behavioral feedback or only consider behavior adaptation by susceptible individuals can 
potentially yield misleading parameter estimates or omit significant features of the disease trajectory. 
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1. Introduction 


In spite of advances in treatment and prevention, sexually transmitted infections (STIs) remain 
endemic worldwide. The CDC estimates that 20 million new cases occur annually in the United States 
alone [1], incurring a total cost of $16 billion for treatment and care. Globally, treatable STIs are 
responsible for approximately 500 million new cases per year [2f, while an estimated 35 million individuals 
currently live with HIV. These statistics underscore the importance of understanding the dynamics that 
drive and sustain STI transmission. To this end, mathematical epidemiology has made substantial 
progress investigating the role of contact patterns such as age-structure, sexual networks, and levels of 
sexual activity However, many open questions remain in understanding the feedback relationship 

between behavioral change and disease dynamics. 

From a behavioral standpoint, sexually transmitted diseases are noteworthy as they require a direct 
and intimate interaction between individuals. As a result, many common preventative measures, such 
as condom use, are not determined unilaterally [5i, 6j. In addition, assuming individuals form preferences 
over protective behaviors based on the associated costs and benefits, we would expect these behaviors 
to respond to the risk of infection as an outbreak progresses [7HS] ■ While changes in risky or protective 
behavior can amplify or suppress outbreaks, adoption of these behaviors can in turn be driven by the 
spread of disease, as demonstrated by increased testing and condom use among men who have sex with 
men in response to the HIV outbreak in the US [10,11]. 

Methods from game theory provide a framework with which to capture this feedback, grounded in 
well-established mathematical and economic theory. The resulting economic- epidemiological models 
can explicitly represent the decision process of individuals either in direct interactions (e.g. sexual 
encounters) or population interactions (e.g. vaccination behavior) [T2H18] . Including the effects of 
behavioral change on STI dynamics has been primarily motivated by the HIV epidemic among men who 
have sex with men (MSM) in the 1980s and 1990s, but this modeling approach is relevant for the study 
of other pathogens and communities as well. Indeed, in 2013 the WHO highlighted the need to study 
behavioral change in order to design effective interventions [2j. 

Many economic-epidemiological models rely on two behavioral assumptions that are worth consid¬ 
eration. The traditional game theoretic framework assumes that all actors are fully rational, responding 
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optimally at every stage of play |T5J. While convenient, the rationality assumption remains a subject of 
debate in economic literature. Empirical studies note circumstances in which individual behavior appears 
to depart from a strict payoff maximization foundation [201422] . In addition, most models (e.g.[23j and 
[24]) assume that only susceptible individuals make choices regarding protective behavior. This assump¬ 
tion is the result of representing the cost-benefit calculus of individuals as a tradeoff between various 
private costs of protective behavior and the risk of infection. However, infected individuals may have 
incentives to reduce contact or use protection as well, motivated by altruism, self interest, or other 
factors [251E6]. This has been the focus of several intervention strategies in practice [271 22]. 

In this paper, we propose a model of combined behavioral and disease transmission dynamics that 
uses the outcome of sexual interactions between susceptible and infected individuals to determine the 
effective contact rate for a mass action model of disease transmission. The combined model bears 
some similarity to the behavior-disease model proposed in JT3j. There are, however, several critical 
distinctions. We use a deterministic ODE framework, while our game-theoretic model collapses the 
protection-use game to a single interaction instead of a multi-stage negotiation. In addition, similar to 
[14] and [T7] we use evolutionary dynamics to represent the process of behavioral change over time, 
allowing for non-optimal but potentially more realistic behaviors. Unlike the inductive reasoning game 
developed by Breban et al. [29], our behavioral dynamics only explicitly consider the current state 
instead of a history of actions. While this approach loses some realistic features, it still allows us to 
relax the assumption of full rationality while also providing a convenient mathematical formulation for 
the combined model 30 [3T]. 


2. Model 

The Susceptible-Infectious-Susceptible (SIS) model has been studied extensively as a simplified 
representation of bacterial sexually transmitted diseases [321134] . The model equations are 


S = jl- psi, 
i = psi- 7 /, 


( 1 ) 
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where S is the fraction of the population which is susceptible, / is the fraction infected, /3 is the 
effective contact rate; the product of the rate of sexual partner acquisition and the probability of disease 
transmission from an infected to a susceptible partner, and 7 is the rate of recovery or treatment. The 
basic reproductive ratio is 

Ho = ^ ( 2 ) 

7 

The disease free equilibrium (DFE) occurs if IZq < 1. Otherwise the endemic prevalence is 


1 1 P' 


(3) 


In order to capture the potential for individuals to adapt their protective behaviors over the course of 


an outbreak, we define a Bayesian game [35] between a pair of players. In Appendix A we give a brief 
overview of the definitions and structure of Bayesian games, with more complete treatments given in 
[G, MM- The payoffs for the game depend on the disease states of both players, which are considered 
private information. Players must infer the type of their partner, reflecting realistic uncertainty about 
serostatus [26] I36H39] , Consistent with the notation for the SIS model, a player may be one of two 
types chosen from the type space 0 = {S,I}. Each player chooses between using protection (P) or 
no protection (U) for a given sexual encounter. If both players select the same action, the outcome 
of the game is the same as the chosen action. We assume that if both players select different actions 
the encounter does not proceed and the effective contact rate for the pair of players is 0. In order 
to characterize the strategy space for this game, it is convenient to use the type-contingent notation 
c Tj(0i ) G {P,U} to denote the action player i would choose if she was of type 0 t . A complete strategy 
for a player then specifies a pair of actions c tj := = S)<jj{Qi = I) G { PP , PU, UP, UU}. Without 

loss of generality, the type-dependent payoff (utility) to player 1 for a given pair of actions and types 
is written u\{aj(Q\), ak(&2), &i, # 2 ) (where the first entry specifies player l's action assuming type 8 1 
and the second entry player 2's action assuming type 6 * 2 ). Then player l's overall expected payoff for a 
particular strategy profile (i.e. for a pair of type-contingent strategies Oj and 07 . for each player) is the 
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Figure 1: The two variants of the protected sex game for the cases a: Oi = 6j and b: 0,; yf Oj. 


double expectation of ui{aj(0{), <7k(02), 9i, O 2 ) over both players potential types |T9], that is, 


E[u 1 (a j (di),ak(d 2 ), 0 1 ,0 2 )\ = ^ Pr(9i)[^2p 1 (6 2 \0 1 )u 1 (crj(d 1 ), a k (0 2 ), Ox, 6 2 )\ (4) 

6»1 02 

where E[u\{aj{6i), a k {0 2 ), Q\, 0 2 )] is sometimes written more simply as E[u\{aj,ak)}- The prior 
probabilities for each player being of either type, Pr{6 £ {S,I}), are given by the distribution of 
susceptible and infected individuals in the population. We assume that partner selection is not assortative 
in disease type, so the belief for player 2 being type 0 2 given player 1 being type 6\ is pi(0 2 \6{) = Pr(0 2 ). 

Fig. a shows the two possible payoff matrices for the protected sex game that determine the type- 
dependent payoff terms Ui(aj(0i), cr k {0 2 ), 0\, 0 2 ) in Equation |4j For a concrete example, we compute 
the expected payoff to player 1 of playing strategy Oj = PU if player 2 picks a k = UU when 30% of 
the population is infected and we choose a = 1, b = 0.75, c = 0.5, d = 0.25. From Equation |4]we have 

E[ Ul (PU, UU)} = 0.7(0.7 x 0.5 -F 0.3 X 0.25) 

+ 0.3(0.7x0.25 + 0.3x1) (5) 

= 0.44. 


In general, we suppose the specific payoff entries satisfy a > b > c > d. Thus, individuals prefer 
unprotected sex to protected sex with a partner of the same disease status, preferring either to action 
pairs resulting in no sexual encounter. However, individuals prefer protected sex to all other outcomes 
when their partner is of a different disease type. This is intended to capture the notion that both 
susceptible and infected individuals have some incentive to avoid infection or transmission respectively. 
We explore the alternate case where infected individuals do not distinguish between susceptible and 


infected partners in Appendix E 
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Consequently, for the protected sex game, the time-varying payoff matrix M is a 4 x 4 square matrix 
with elements given by Eq. Q. We can compute the Bayes-Nash equilibrium for any given choice of 
a,b,c,d, S, I. One case of particular interest is at the disease free equilibrium, I* = 0. Here the 
game reduces to the 2x2 game with payoffs as in Figure [Ta. This game has two pure strategy Nash 
equilibria, (P. P) and (U,U). However, the (U,U) equilibrium in the reduced game is both payoff and 
risk dominant. Since the a(0i = /) actions do not contribute to the expected payoff, any mixture of 
the type-contingent strategies UP and UU (similarly PP and PU ) is a Nash equilibrium in the full 
model. We will revisit this point shortly. 

Since we are ultimately interested in the dynamics of behavior at more than a single disease state, 
we use methods from evolutionary game theory to couple the dynamics of strategy change over time 
to the disease trajectory. This allows us to close the feedback loop between behavior and disease 
dynamics. In particular, we use replicator-mutator dynamics with a linear fitness function [311 40J. 
The distribution of strategies in the population using the replicator-mutator approach depends both on 
the existing distribution of strategies and their expected payoffs, as well as a small degree of random 
strategy choice. This allows us to capture the notion that individuals may not respond immediately or 
strictly optimally to the presence of an outbreak, as payoff-suboptimal strategies may remain frequent 
in the population for some time. In addition, the possibility of mutation prevents any strategy from 
becoming extinct. In the behavioral context, this can be thought of as allowing for a small amount of 
random choice. This feature is particularly important as the replicator equation has stable steady states 
at strict pure strategy Nash equilibria, such that behavior cannot change after fixation on a particular 
Nash equilibrium even if disease conditions change. Both players in the protected sex game have the 
same action and type sets, so the relative frequency of the j th strategy in a large population is f(<jj,t). 
The replicator-mutator equation for a given strategy aj is 

f{°j) = f((T k ,i)qkj^((Tk,t) - 4>(t)f(<jj,t)}, (6) 

k 

where q^j is the probability that an individual playing aj switches to <r fc , <f)(aj,t ) = E[u(<Tj,aj c )]f(<7k,t) 

is the fitness of aj , and s G [0,oo) is a scaling term that determines the 
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speed of behavior change. In vector-matrix form, this can be written as 


f = s[D f Q(Mf) - (f T Mf)f], (7) 

where M is the 4x4 payoff matrix with elements rrijk = E[ui(<jj, <7*)], f = ( f(PP , t), f(PU, t), f (UP, t ), 
f(UU,t)), Q is the mutation matrix and Df = diag( f). Note that when Q = I, the evolutionary dy¬ 
namics are equivalent to standard replicator dynamics. However, for our simulations we use a mutation 
probability // = 0.03 so 

« = (i-|rtr + f-i- (8) 

In a large population, the aggregate effective contact rate is determined by the average outcome over all 
pairs, so (3 = where fib is the baseline effective contact rate, Sjj = S x (f(UP,t) + f(UU,t)) is 

the proportion of susceptible individuals playing U, and Ijj = I x (f(PU,t) + f(UU,t)) is the proportion 
of infected individuals playing U. 

The combined model equations with evolutionary behavioral dynamics can be written as 

S = 7 / - (3bSuIu, 

i = PbSulu - 7/, (9) 

f = s[D f Q(Mf) - (f r Mf)f], 

The system above has six compartments. Of these, only four are strictly necessary as S + 1 = 1 and 
= 1- In a completely susceptible population, action pairs where both susceptibles play P 
or U are Bayes-Nash equilibria. In order to restrict the domain of possible initial conditions, we take 
the risk dominant equilibria, or mixtures between UP and UU. This can be interpreted as individuals 
in a disease free state preferring unprotected sex because it best hedges against uncertainty in partner 
actions. While the a(Oi = I) actions are never realized at the DFE, the type-contingent strategy 
framework implies that individuals must be able to specify an action that they would take if they 
became infected. The effective contact rate for a single infected individual is (f(UU,0) + f(PU,0))(3b, 
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where f(UU,0),f(PU,0) G [0,1]. As a result, the basic reproductive rate of the combined model is 


= (/(C/f/,0) + /(PP,0))& 

7 

Thus, IZo depends on the behavioral initial condition, unlike the standard SIS model. The value of the 
above expression will lie within the interval [0, (3bh\- The payoff and risk dominance of the (U, U) Nash 
equilibrium at the initial disease-free steady state (discussed above) suggests that /(PC/,0) is likely to 
be low, and indeed for the remainder of this paper we will assume that the initial infected plays the 
strategy UU, while the underlying susceptible population at t = 0 will be assumed to take on a mixture 
of the UU and UP strategies. It is nonetheless of interest to note that the dependence of TZq on the 
initial behavioral conditions suggests the potential for both high TZ o with disease extinction and low TZq 
with endemic prevalence: if /(PC/,0) ~ 1 it may be possible to have high TZq but still have disease 
extinction due to all susceptibles playing the protected strategy, while conversely if TZq is low due to 
low frequency of UU and PU but otherwise has a high contact rate, it may be possible to generate 
endemic prevalence even for TZq < 1 (as explored further below). 

3. Results 

Due to the addition of the replicator-mutator equations, it is not straightforward to solve for the 
steady states of the combined model analytically. Instead, we conducted a range of numerical simulations 
to investigate the effects of behavioral dynamics using Python 2.7 with Numpy 1.9, Scipy 0.15.0, and 
Matplotlib 1.4.2. In order to explore a range of outcomes generated by the combined model, we focused 
on both long-term and short-term dynamics. Given the addition of behavior change, the short term 
dynamics provided useful mechanistic insights to explain the steady state properties of the model. Unless 
noted otherwise we used a = 1,6 = 0.6, c = 0.4, d = 0.2 for the values of the outcome payoffs. In 
general the qualitative features of the model do not change substantially for different payoff values 
provided the overall scale and ordering are preserved. Unless otherwise specified, the disease initial 
conditions were 5(0) = 0.99, 1(0) = 0.01. In the first section below, we consider fixed initial behavioral 
conditions at f(UU, 0) = f(UP, 0) = 0.5, while in the following section we consider the effects of 
varying the initial behavioral conditions (for f(UU, 0) between 0 and 1). 
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Figure 2: Simulated trajectories for the combined model with Rq = 1.2. Disease parameters /3b = 1-2,7 = 0-5 
with behavioral parameters a = 1,6 = 0.6, c = 0.4, d = 0.2. (a) Disease prevalence, (b) Type-contingent 
strategies, (c) Fractional contact reduction, (d) Strategy fitness. 
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3.1. Model Dynamics 

We began by inspecting an example set of model trajectories, shown in Fig. [2] (with 06 = 1.2,7 = 
0.5, and /(C/C/,0) = 0.5 to give IZq = 1.2), to illustrate the interactions between behavioral changes 
and disease dynamics. Comparing the prevalence and contact reduction trajectories demonstrated the 
impact of the adaptive behavioral dynamics on the progression of the simulated outbreak. In this 
example the outbreak was relatively small, so susceptible individuals did not have a large incentive to 
use protection. However, the initial increase in disease prevalence favored strategies where infected 
individuals used protection, which in turn decreased the force of infection enough to halt the initial 
outbreak. However, as prevalence decreased, unprotected strategies became less costly and contact 
reduction decreased again. Consequently there was a small secondary outbreak. The timing of the 
behavioral response also appeared to be influential. While contact reduction tracked prevalence over 
time, it did so at a delay. Thus the level of contact reduction was still relatively high at the onset of 
the second outbreak, preventing a large secondary peak. This preliminary exploration suggested that 
both the infectiousness of the disease and the speed of behavioral adaptation play important roles in 
the overall model dynamics. 

To explore this issue further, we defined a timescale parameter s for the behavioral dynamics, given 
as a scaling factor on the payoff matrix A in Eq. 0 From a behavioral perspective, this parameter 
can be thought of as controlling the speed of adaptation in the population. We then evaluated the 
long-term dynamics of the model in steady state for a range of values of IZq and s, shown in Fig. [3] 
(with all remaining parameters as in Fig. [2]). In the triangular region about IZq = 2.7 and the canyon 
between 3 < IZq < 3.5 the damped oscillations became sustained with a substantially higher average 
steady state prevalence. The triangular region of oscillations gives way to another endemic steady state 
region as IZq increases, followed by the thin band of oscillations and eventually disease extinction in the 
triangular region along the right edge. 

Figure [4] and Supplementary Figure [B~8] provide a more detailed examination of the dynamics in 
Fig. [3] as IZq increases for s = 1. For IZq sufficiently above 2, the early increase in contact reduction 
by infected individuals was not sufficient to stop the spread of the initial outbreak (JjJj). As a result 


prevalence increased, causing infecteds to begin to switch back to unprotected strategies (B.8 b). 
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Figure 3: -D- Long term dynamics of the combined model with increasing values of and the behavioral scale 
parameter s. IZo was adjusted by varying the parameter /3b in Eq. ( flO| ). (a) The amplitude of steady state 
prevalence oscillations, (b) The average prevalence at steady state. 


However, at higher prevalence susceptible individuals had an incentive to use protection, compensating 
for the behavior of the infected population. As in the initial example, this overall increase in contact 
reduction reduced the force of infection and halted the growth of the outbreak. The magnitude of this 
response appeared to determine the long term dynamics of the model. When prevalence was reduced to 
a moderate level neither protective nor unprotected strategies were able to gain a stable foothold in the 
population, leading to sustained oscillations in both contact reduction and prevalence (pffr, |B.8 b). A 
stronger contact reduction response pushed the disease to a level where infected individuals again had 
a strong incentive to use protection. The fitness advantage of this strategy was sufficient to survive the 
second outbreak, leading to an endemic steady state ([4]:, |B.8| c). Notably, when TZq was high enough 
that the resulting behavioral response nearly drove the disease to extinction, sustained oscillations were 
again possible (|4|j, |B.8 d). Subsequent outbreaks in this scenario were contained by increasing contact 
reduction from a higher baseline due to the initial outbreak, but grew rapidly due to the infectiousness of 
the disease. Finally, when the initial outbreak was extremely large, essentially all remaining susceptibles 
(including those who recently recovered) played protective strategies, driving the disease to extinction 
?8] e). Since the boundaries of both oscillatory regions shifted with both s and IZq, it appeared 
that the timescales of infection and adaptation must align in order to produce sustained oscillations 
(e.g. Figure [4jo and B.8a). In particular, oscillatory dynamics appear only to be possible for relatively 
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Figure 4: Prevalence and contact reduction from the combined model with 7 = 0.5, f(0) = (0.0,0.0,0.5,0.5) for 
all simulations and (a) /3b = TZo = 1.5, (b) 1Z 0 = 2.7, (c) IZq = 2.9, (d) 1Z 0 = 3.3, and (e) 1Z Q = 3.5. 
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Figure 5: Average steady state prevalence of the combined model for increasing with 7 = 0.5, f(0) = 
(0.0, 0.0,0.5, 0.5). The transition between the disease free and endemic equilibrium occurs at 1 < Tvlo < 1.05. 

slow adaptation speeds, suggesting that the lag induced in the behavioral dynamics contributes to the 
potential for oscillations. 

Using the same parameters as in Fig. [4j in Figure [5] we also examined steady state prevalence near 
TZq = 1, to evaluate whether behavior changes may affect the threshold for generating an outbreak. We 
found that steady state prevalence remained zero even when TZq was greater than 1 , suggesting that 
behavioral dynamics were able to extinguish the disease even when the initial growth rate would have 
generated and epidemic. 


3.2. Behavioral Initial Condition Analysis 

For our next set of simulations, we tested the effect of varying the initial behavioral conditions. 
Figure^ shows the average steady state prevalence across /(C/C/,0) £ [0,1] when /7 = 3 (where 
the remaining population was assumed to take strategy UP as noted above). As suggested by the form 
of n 0t an endemic steady state persisted even for TZo < 1. This counterintuitive phenomenon occurred 


because, while the disease initially declined, it did not immediately become extinct (Fig. C.9). The 
declining prevalence led to increased adoption of C/C/, driving the effective reproductive rate back above 
one. Figure [ 6 |o depicts the effect of different levels of f(UU, 0) on the steady state behavior of the 
combined model for increasing flb/'y. A wide range of dynamics were observed as {3b /7 varied, including 
a single stable equilibrium, bistability, and oscillations. When /(C/C/,0) was low (green line in Fig. [ 6 }o) , 


the model only exhibited damped oscillations, and did not produce the same extinction behavior at high 
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Figure 6 : The effect of changing the initial condition /(C/C/,0) with 7 = 0.5. (a) Steady state prevalence and 
7?-o for /3b/'y = 3 (Dashed line indicates where 7\t 0 = 1). (b) Steady state prevalence of the combined model as 
Pb/'y increases for three choices of /(C/C/,0). In regions of (3b /7 where sustained oscillations occurred, the solid 
black line shows the value of the crest while the dashed line shows the value of the trough. 


j3b /7 as in the case where /(C/C/,0) was substantially larger than zero. For /(C/C/,0) = 0.5 (blue dotted 
line), the model exhibited the same bifurcation pattern as in the previous section. At high /(C/C/,0) 
(black line) the behavior was similar, however the endemic prevalence before the first oscillatory region 
was substantially greater. From this, /(C/C/,0) appeared to act as a switch between possible steady 
state regimes of long term dynamics. 

For the constant steady state regions, the endemic prevalence took one of two values for a given 
Pb/l depending on /(C/C/,0), similar to Fig. [ 6 ^. More broadly, for any fixed /3b/ 7 , there were two 
basins of attraction corresponding to either the /(C/C/, 0) = 0 or /(C/C/, 0) = 1 case, where the unstable 
equilibrium dividing the two regions depended on the value of /3b/ 7 - This is illustrated by the /(C/C/, 0) = 
0.5 curve in Fig. §3. which switches between existing in the /(C/C/,0) = 0 and /(C/C/,0) = 1 basins 
around /J &/7 = 5. There are also regions of Fig. [ 6 }) where only a single steady state exists regardless of 
the value of /(C/C/,0), at very low values of /3b /7 (left portion just above (3b /7 = 2) and larger values 
(between the oscillatory regions). 

As noted in the previous section, there were also multiple regions where the long term dynamics 
showed disease extinction even though IZq > 1. For /(C/C/,0) = 0.5 and 1, the left corner of Fig. & 
shows extinction with TZq > 1 , and the right corner shows the same extinction discussed in the previous 
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section for the /(C/C/, 0 ) = 1 basin at high values of Pb/l- 


3.3. Comparison with Alternate Models 

Finally, we considered how the dynamics of the model compare to other potential models of the 
disease dynamics, to examine how neglecting the behavioral dynamics may alter model forecasts of the 
epidemic trajectory or affect estimates of key epidemiological parameters such as TZo. As a first example, 
we chose a fixed contact rate SIS model parameterized such that the IZq of both models was identical. 
While the trajectories were similar at very early times (at which the epidemic growth rate can still be 
characterized by IZq), the combined model quickly diverged due to the reduced effective contact rate, 
and equilibrated at a substantially lower endemic level (FigurejT^). Similarly, we computed IZq = 1 _ 1 J 
naively from the steady state prevalence of the combined model in Figure [7^, without accounting for 
a time-varying contact rate. The estimated value of 1.03 was substantially lower than the true IZq of 
1.2 for the combined model. This difference can be even more severe when considering endemic steady 
states at higher IZq's (such as the area to the right of the triangular region in Figure [3]). For example, 
the endemic steady state in Figure [4J: yielded an apparent TZq of 1.21, when the underlying IZq for the 
combined model was 2.9. Given the simplicity of the SIS model, these discrepancies were not surprising. 
Nonetheless, they highlight some pitfalls of neglecting the effects of behavioral dynamics. 

However, as the SIS model is known not to be able to produce oscillations, it might not be a likely 
choice given data that came from real-world infection dynamics similar to the combined model. Due to 
the delay between recovering from infection and returning to susceptibility, the SIRS model is capable 
of producing damped oscillations, and so might make a reasonable initial guess for the model structure 
if the length of immunity was unknown. To evaluate how a more realistic but misspecified model might 


affect parameter estimation, we fit the SIRS model (Appendix D) to the prevalence trajectory of the 
combined model with TZo = 1.2, using least squares with Nelder-Mead optimization, as shown in Figure 
[7}x While the best fit SIRS model conformed well to the target trajectory, the resulting parameter 
estimates included both a four fold decrease in the baseline contact rate and a substantial increase in 
the average waiting time until a recovered individual becomes susceptible again. The TZq for the best fit 
SIRS model was 1.35, a 12.3% increase compared to the behavioral model. In this case, the feedback 
between behavior and transmission may make interventions more effective than would appear to be the 
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Figure 7: Comparison of simulations of the combined model other compartmental models, (a) SIS model contact 
rate (3 = 0 . 6,7 = 0-5 to match the R 0 = 1.2 of the combined model, (b) Prevalence trajectories for the 
combined model and an SIRS model fit to the combined model prevalence simulated with TZq = 1.2. SIRS 
best fit parameters: (3 = 0.274,7 = 0.203,(5 = 0.022. (c) A variant combined model in which only susceptible 
individuals adapt their behavior fit to the simulated trajectory from the full model with /3b = 1 . 8,7 = 0-5 and 
TZq = 2.7. The best fit parameter (3£ = 2.84 for the reduced model, giving an TZq = 5.68. 
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case if we had assumed a fixed contact rate. 


Finally, most previous work combining game theory with transmission models to model contact 
reduction has focused on the adaptive behavior of susceptibles [T4j [18, 231 [24]. This is a natural 
formulation for vaccination, but may be less suitable for more general contact reduction behaviors. 
Thus, to test the effect of modeling adaptive behavior by both susceptible and infected individuals, we 
used a reduced variant of our combined model in which infected individuals always select the action U. 
As with the SIRS model, we fit this reduced model to a simulated trajectory from the full model. In 
this case, the reduced model overshoots the initial outbreak curve, but equilibrates to a similar endemic 
equilibrium to the full model (Figure]?]:). However, the best fit parameter and IZq for the reduced model 
were nearly two times higher than the true values for the full model. For a wider range of parameter 
values of the full model, the reduced model generally failed to provide qualitatively accurate fits. In 
addition, the reduced model did not produce the same range of dynamics as the full model. We give a 


full description of the reduced model in Appendix D as well as an expanded set of fits and simulations. 


4. Discussion 

To explore the feedback between behavior and disease dynamics, we developed a model of disease 
transmission with adaptive behavior among both susceptible and infected individuals. Numerical simu¬ 
lations illustrated the effect of behavior-disease feedback on model dynamics and inferences about STI 
transmission (Figure]?]). We found that phase transitions between damped and sustained oscillations 
occurred at an intermediate transmission rate and again at a higher transmission rate (Figure [3]). In¬ 
creasing the adaptation rate parameter s reduced the range of transmission rates at which sustained 
oscillations could occur. 

Disease extinction occurred both for low and high transmission rates (Figures [3] [6]). At lower 
transmission rates, adoption of protection by infected individuals lead to extinction above the typical 
TZq = 1 threshold (Figure]?]). This suggests that behavior-disease feedback can create a herd-protection¬ 
like effect if the infected population can be reduced sufficiently to incentivize the use of protection among 
that group. At high transmission rates and low-to-mid range levels of initial risky behavior (/(C/C/,0)) 
in the infected population, we observed the counterintuitive result that a large initial increase in risky 
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behavior (C/C/ and PU strategies) among infected individuals preceded disease extinction. However, 
because susceptible individuals were also switching to protective strategies (PP and PU) rapidly, the 
combined effect of both susceptibles and infected adopting PU resulted in sufficient growth of PU to 
drive the disease to extinction as the effective number of susceptible individuals remained small and 
all were largely playing protective strategies (Figure [4 B.8). In fact, both the high TZq and low TZq 
extinction phenomena (for TZq > 1 but still low) appear to be driven by a reduction in the effective 
size of the ‘bottleneck’ population group. The low TZq extinction (Figure [5]) seems largely due to 
increasing infected protective strategies, while the high RO extinction seen on the right side of Figure[4] 
is largely due to increasing susceptible protective strategies. In each case, the low/high prevalence makes 
infecteds/susceptibles a bottleneck group, where sufficient protective response in that group appears to 
result in disease elimination. The feedback between disease and behavior thus often constrains endemic 
outcomes to scenarios where a disease is only moderately infectious, for a range of realistic levels of 
initial risky behavior (Figure [6^). 

However, simulations predicted a more pessimistic outcome for changes in the initial strategic 
distribution. In general, the initial level at which infected individuals used protection did not determine 
whether the disease would become extinct. As long as the value of ^ 6/7 was sufficiently large, any 
/(C/C/,0) led to an outbreak and endemic disease (Figure [6]). Instead, the initial distribution acted 
as a switch between two steady state regimes: one in which initial risky behavior by infecteds (high 
/(C/C/,0)) lead to sustained oscillations and higher endemic prevalence could occur for values of q, 
below the first oscillatory region, and another in which initial protective behavior by infecteds (low 
/(C/C/,0)) yielded damped oscillations and generally lower endemic prevalence. The exception to this 
pattern is the region where fib/'y is large. As noted above, in this region extinction occurs when early 
risky behavior by infected individuals is then countered by widespread use of protection by susceptible 
individuals. However, early widespread protection use by infecteds leads to a smaller outbreak that 
is slowed, but not eliminated by susceptible protection use. This complicates the interpretation of 
TZq for the combined model, as the standard endemic threshold property often does not apply. For 
interventions it is thus important to understand the interaction between early behavioral patterns, the 
transmission rate, and treatment levels. Interventions based on these factors may yield counterintuitive 
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results depending on the other factors, making comprehensive data collection critical when optimizing 
intervention strategies. 

It is also interesting to note that while an analysis of the Bayes-Nash equilibria of the protected 
sex game at DFE indicates that the susceptible-unprotected strategies (U,U) are risk dominant, our 
combined model suggests that the actual initial distribution will depend on a population’s previous 
experience with disease. That is, when the outbreak is large (i.e. for high TZq, Figure [4]), adoption 
of protective strategies results in disease elimination. However, the increased protective strategies by 
susceptibles mean that the model switches to the attraction basin of the non-dominant (P,P) Nash 


equilibrium as it goes to the new disease-free steady state B.8 (or as close as possible given the mutation 
rates in the replicator-mutator equation). By contrast, below the threshold for high Ro extinction, the 
model tends to reach a steady state close to the susceptible-unprotected (U,U) Nash equilibrium 
since endemic prevalence is low enough to favor susceptible-susceptible contact. This suggests that in 
general, the disease-free strategy balance between protected and unprotected depends on the magnitude 
of previous outbreaks, with larger outbreaks yielding disease-free protective behaviors even though they 
are neither payoff nor risk dominant. As a result, post-outbreak surveillance and control are crucial to 
prevent recurrent outbreaks due to reintroduction. 

Based on comparisons between the combined model and similar disease models with simplified 
contact dynamics, predictions from fixed-contact rate models may omit important dynamical features or 
yield misleading parameter estimates (Figure]?]). Similarly, if an STI model is misspecified, the common 
method of computing IZo from endemic prevalence when a disease is assumed to be at equilibrium 
[4lj is not appropriate and yields inaccurate estimates. One advantage of using a game theoretic 
framework to model interactive behavior is the flexibility of ordinal utility. Since protective measures 
can be used by susceptible or infected individuals to prevent infection or transmission, respectively, we 
used a preference order that gave individuals in both disease states an incentive to use protection if 
paired with a partner of the opposite state. Our simulation results suggest that disease dynamics differ 
in this scenario as opposed to when only susceptible individuals adapt their behavior, and indeed model 
fits using the reduced model considering only susceptible behavior change often yielded both incorrect 


TZq estimates and dynamic trajectories (Figure D.ll). The role of infected behavior dynamics was 
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particularly notable in the early stage of outbreaks when disease was not prevalent enough to induce 
susceptible individuals to use protection. In addition, the shift between infected and susceptible adoption 
of protective strategies appeared to drive sustained oscillations. To more thoroughly consider the effects 
of infected behavior dynamics, we also considered an alternative preference structure in which infected 
individuals no longer explicitly prefer protected sex with susceptible individuals, but rather have the 


same preferences regardless of their partner’s disease status (Appendix E). This model only exhibits 


one oscillatory region and a higher endemic prevalence than the original model (Fig. E.12). However, 
extinction still occurs at high IZq, confirming that susceptible protective strategies drive the high-7£o 
extinction phenomenon. 

In order to focus on the effect of adaptive protective behavior, our combined model simplifies many 
other realistic factors that contribute to STI transmission. The protected sex game does not include 
explicit negotiation, which likely biases the effective contact rate downward, since negotiation or coercion 
could lead to unprotected sex even when the initial action pair is (U. P). The combined model presented 
here still uses mass action assumptions to determine interactions between individuals. Although more 
complex contact patterns are known to influence transmission dynamics, we model a homogeneous 
population to focus on the effect of time-varying behavior. While this may not be completely implausible 
in a highly active group such as MSM frequenting bathhouses, it is almost certainly a poor representation 
of the manner in which individuals form sexual partnerships. Typically, adding contact heterogeneity 
increases the TZq of a model, so we might expect more rapid early outbreak growth for a wider range 
of parameters in our model, potentially expanding the regions in which sustained oscillations or even 
disease extinction occur. Our general framework, however, is amenable to extensions to expand the 
state space of the disease model to represent more complex natural histories or population structures 
as with standard transmission models. In addition, it is possible to model the combined process on a 
contact network or using a stochastic framework. One particular extension that could yield insight into 
spatial patterns and the effect of local information would be a simulation of the model on a regular 
lattice, where individuals could use either a global or local prior to estimate the probability of their 
partner's type. 

A particularly compelling consequence of developments in economic-epidemiological models is the 
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potential to estimate more complex behavioral parameters using traditional surveillance data. In the 
context of our combined model, this would result in estimates for the relative payoff values ( a,b,c,d ), 
which capture useful information about the perceptions of populations facing STI outbreaks. While 
subject to as yet unknown identifiability properties, this manner of parameter estimation could provide a 
valuable link between game theoretic methods and the extensive empirical literature on the epidemiology 
of STIs. 
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Appendix A. Bayesian Games 

This section provides a brief overview of static Bayesian games of incomplete information. Those 
desiring a more complete treatment may refer to |19j.|42|.and [35]. 


Appendix A.l. Definition 

A normal form n-player symmetric static Bayesian game includes 

• A set of players N = {1, 2,..., n}. 

• Actions a* G A for each player. 

• Types Qi e 0 for each player. 

• Pure type-contingent strategies a : 0 —> A. By convention, we use the notation 

= a A 9 \) G A 9 !)- .Oj( 9 ™) to represent the jth strategy for a finite type space. 

• Belief distributions p t where pi{9-i\9i) is the conditional distribution on the types of other players 
given player % s type, —i denotes the set of players except i. 

• Type-dependent payoffs Ui : A n x 0 n —> M for each player. The expected type-dependent payoff 
(E[ui(ai,...,a n \9i)]) gives the average payoff over player i’s belief regarding the other players’ 
types conditional on player i’s own type. E(ui{a\,a n , 9\,9 n )) is player i’s unconditional 
average payoff over all n player types, sometimes denoted E(ui(aj 1 , crj n )) in terms of strategies. 

In general, a profile of types, actions, or strategies is defined to be a listing of the particular 
types/actions/strategies (respectively) assigned to each player. While the actions in a static Bayesian 
game take place simultaneously, it is useful to break the game down into stages as follows 

1. Nature chooses a profile of types (9\, 02, ■ ■., 9 n ) from the common prior distribution. 

2. Each player picks a strategy (jj i 

3. Each player learns only his type A;. 

4. Using Bayes’ theorem and the common prior, each player forms beliefs pi{0_i\0i) over other 
players' types. 
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5. Players choose actions simultaneously according to their strategy (ai = <7j, (0\ )) to form a profile 

(«1 1 j • • • j On)- 

6 . Players receive their payoffs «j(ai, a n , 9±, #„) based on the action profile and the type profile. 


Appendix A.2. Bayesian Nash Equilibrium 

The strategy profile (cr^,..., cr* J is a Bayesian Nash equilibrium if for all players and all types, a* 


ji 


satisfies 


£ 

0-t 


Pi(6-i\0i)ui(a*.(0i),<j*_.(6-i),9i,theta-i) > ^2pi(9-i\9i)ui(a ki (9i), cr*_.{9-i), 9 t , 9-i) (A.l) 


0-i 


for any a k . / cr*.. Equivalently 


0j 0-i 9i 0_i 

(A.2) 
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Appendix B. Behavioral Trajectories for Simulations in Figure |4] 
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(a) (b) 




(c) (d) 



(e) 

Figure B.8: The distribution of type-contingent strategies as a function of time with f(UU,0) = 0.5,7 = 0-5 as 
in Figure [4] 
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Appendix C. Endemic Prevalence at TZq < 1 



Figure C.9: Prevalence and contact reduction from the combined model for /3b = 2.1 ,7 = 0.5, f{UU , 0) = 0. 


Appendix D. Alternate Models 

Appendix D.l. SIRS Model 


S = 5R — PSI 

i = psi - 7 / (D.l) 

R = jl- 5R 


Appendix D.2. Susceptible-only Behavior Change 
For the reduced model, we use the fitness functions 

<KPs,t) = Pr(S)[f(P s ,t)u 1 (Ps,Ps,S,S) + f(Us,t)u 1 (Ps,U s ,S,S)\ 
+Pr(I)[f(P I )u 1 (P s , Pi, S, I) + /(U^u^Ps, Uj, S, /)] 

(D.2) 

<KUs,t) = Pr(S)[f(Ps,t) Ul (Us,Ps,S,S) + f (Us, t) Ul (Us, Us, S,S)} 
+Pr(I)[f(P I )u 1 (U s , Pi, S, I) + f(Ui) Ul (Us, Ui, S, /)] 
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where the notation P$ denotes a susceptible player choosing P (similarly for C/ 5 , Pj, Uj) and f(Uj) and 
F(Pi ) are fixed over time. Since only susceptible players change their strategy, f(t) = (f{Ps, t),f(Us, t)), 
the mutation matrix for this game is 2 x 2 and we can use the two-dimensional system 


i = -f%Suf{ui)i + 7 / 

f(Us) = qu s p s 4>(Ps,t)f{Ps,t ) + qu s ,u s ( l ) {Us,t)f{Us,t ) - <f>f(Us,t ) 


(D.3) 



model compared to the full model. 



Figure D.10: Prevalence dynamics for the reduced model from Section [2] where only susceptible individuals adapt 
their behavior. For the simulations above 7 = 0.5, f(0) = (0.0,1.0). 


Appendix E. Alternate preferences 

We examined a case where infected individuals prefer unprotected sex over protected sex regardless 
of partner type, a phenomenon that has been observed empirically [43J■ For susceptible individuals, the 
type-dependent payoffs from the alternate game are the same as in Figure[l] Without loss of generality, 
for an infected-type player 1 paired with a susceptible-type player 2 , the type-dependent payoff matrix 
is 
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Player 2 

P U 


Player 1 


P 

U 


b 

c 

c 

a 


Figure [El2] shows the long-term behavior of the alternate model for a range of TZq and s. This model 
only exhibits one oscillatory region and a higher endemic prevalence than the original model. However, 
extinction still occurs at high TIq. 
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(a) (b) 




(c) (d) 




(e) (f) 

Figure D.ll: The reduced model of susceptible-only behavior change fit to simulated trajectories from the full 
model with 7 = 0.5, (a) Full model (3b = 1.2, best fit (3 J = 0.542, (b) full model (3b = 1.8, best fit (3(( = 3.41, 
(c) full model (3b = 2.4, best fit /3j* = 2.84, (d) full model (3b = 2.7, best fit (3^ = 1.99, (e) full model (3b = 3, 
best fit (31 = 2.36, (f) full model (3b = 3.6, best fit (3(( = 2.82. 
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Figure E.12: Long term dynamics of the combined model where infected individuals always prefer unprotected sex 
for increasing values of and the behavioral scale parameter s. (a) The amplitude of steady state prevalence 
oscillations, (b) The average prevalence at steady state. 
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